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ABSTRACT 

Spiral arms that emerge from the ends of a galactic bar are important in interpreting ob¬ 
servations of our and external galaxies. It is therefore important to understand the physical 
mechanism that causes them. We find that these spiral arms can be understood as kinematic 
density waves generated by librations around underlying ballistic closed orbits. This is even 
true in the case of a strong bar, provided the librations are around the appropriate closed or¬ 
bits and not around the circular orbits that form the basis of the epicycle approximation. An 
important consequence is that it is a potential’s orbital structure that determines whether a bar 
should be classified as weak or strong, and not crude estimates of the potential’s deviation 
from axisymmetry. 
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1 INTRODUCTION 


Many grand design barred spirals exhibit spiral arms starting at and 
extending out from the end of the bar. How are these spiral arms 
generated? The natural interpretation is that they are driven by the 
bar. The spiral arms and the bar must also be rotating at the same 
pattern speed if the starting points of the arms always coincide with 
the bar’s ends. Such spiral arms are relevant for the interpretation of 
observations of both our Galaxy and external galaxies. For exam¬ 
ple, it is likely that some features present in spectral-line emission 
of atomic and molecular gas in the inner Galaxy, for instance the 
“3kpc arm”, are produced by such arms ( |Dame & Thaddeus|2008| 
e.g.). It is therefore important to have a physical understanding of 
why such spiral arms arise. 

On the theoretical side, the possibility of bar-driven spiral 
arms is now well established. The first investigations used ana¬ 
lytical methods and were carried out by [Feldman & Lin| ( |1973| ) 
and |Lin & Lau|(T975| ). Then, in the following years, [Sanders &| 
[Huntleyl ( |1976| l; Sanders] \ 1911 \ and [Huntley et al.| |l978[ ) per^ 
formed pioneering numerical experiments and discussed the phys¬ 
ical mechanism responsible for the generation of the spiral arms 
in terms of closed orbits. Their work was later extended indepen¬ 
dently by |Wada|([1994[) and|Lindblad & Lindbl^ ( [1994[ ), and more 
recently by Pinol-Ferrer et al.| ( [2012[ ). These authors constructed 
phenomenological analytical models under the epicyclic approxi¬ 
mation with the purpose of understanding the spiral arms. In their 
models gas parcels follow weakly oval orbits, whose major axes 
orientation changes with radius giving rise to kinematic density 
waves a-la Lindblad, and hence to spiral arms. 

An alternative viewpoint which does not consider the spiral 
arms as density waves has also been discussed in a series of pa- 
pers by [Romero-Gomez et ^ < |2006[ [2007| l; [Athanassoula et af] 
( [2009a|b| [2010j l. Their theory, which is more directly applicable 
to stars than to gas, is based on the observation that orbits in the 


vicinity of unstable Lagrangian points can be trapped into invariant 
manifolds whose morphology can reproduce the spiral arms. 

Many authors have noted the presence of bar-driven spiral 
arms in simulations, and discussed them in a more or less descrip¬ 
tive way (s ee for example [Athanassoula|1992[[Englmaier & Ger- 
hard|199^[Bissantz et al.|2003||Rodriguez-Fernandez & Combes 

2008] ). For a recent review see also Section 2.3 of [Dobbs & Baba 

( [20T4| t. Note that in all works cited above the gas is driven by an ex¬ 
ternal potential generated by the bar and is not self-gravitating, and 
therefore the spiral arms are not spiral density waves in the sense 
of [Lin & Shu[ ( [1964] ). We do not discuss self-gravity in this work. 

The viewpoint that spiral arms can be understood as density 
waves is nowadays often assumed, but has never been tested with 
sufficient detail. On modern computers, it is very cheap to run simu¬ 
lations able to test the predictions and the limits of the phenomeno¬ 
logical models cited above. Moreover, the literature only addresses 
the weak bar regime under the epicyclic approximation, and does 
not discuss how the picture should be extended to the strong bar 
case. 


In this paper we investigate the physical mechanism respon¬ 
sible for the generation of the spiral arms and in particular we re¬ 
sume the discussion of how they can be understood as kinematic 
density waves. The structure is as follows. In Sec. we briefly re¬ 
view previous work aimed at understanding in a phenomenological 
fashion the spiral arms in the epicyclic approximation. Then we run 
grid based, isothermal, non-self gravitating 2D hydrodynamic sim¬ 
ulations in an externally imposed rigidly rotating barred potential, 
addressing both the weak and the strong bar case. The numerical 
methods employed are explained in Section In Sec. we com¬ 
pare the results of the simulations with the phenomenological mod¬ 
els available in the literature. In Sec. [5] we discuss in more detail the 
weak bar case, and in Sec.[^how the results for the weak bar case 
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2 Sormani, Binney & Magorrian 


should be extended to the strong bar case. Finally, we summarise 
our findings in Sect.[^ 


2 REVIEW OF PREVIOUS WORK 

In this section we briefiy review previous work aimed at under¬ 
standing the bar driven spiral arms in the epicycle approximation. 
In Sec. |2.I| we discuss the equations describing ballistic closed or¬ 
bits in the epicycle approximation (see for example Sect. 3.3.3 in 
[Binney & Tremaine]|2008| l. In Sec. |2.2| we discuss how various 
authors have modified the ballistic equations in order to describe 
gaseous parcels, and the phenomenological models based on these 
equations. These phenomenological models will be in later sections 
compared with the results of hydro simulations. 


2.1 Ballistic closed orbits in the epicycle approximation 

Consider a rigidly rotating external potential of the form 

( 1 ) 

where Oq is an axisymmetric potential and Oi is a small but other¬ 
wise arbitrary perturbation. The potential is assumed to be rigidly 
rotating at pattern speed flp. 0 are polar coordinates in the rotat¬ 
ing frame, with 0 = 0 corresponding to the positive horizontal axis, 
and 0 increasing clockwise. The equation of motion for a ballistic 
particle in this potential is 

X = — VO + flpX — 2flp (e^ X x) , (2) 

where is the unit vector perpendicular to the plane. The first term 
on the RHS represents gravitational forces, the second term is the 
centrifugal force and the third is the Coriolis force. 

We want to find closed orbits in the above potential under the 
epicycle approximation. In the absence of the perturbation Oi the 
potential is axisymmetric and the only stable closed orbits are cir¬ 
cular orbits. These circular orbits will be described in polar form 
as 


R{t)=Ro, (3) 

9(j) = eo(0 = iif^ (4) 

where 

flf = - flp , (5) 

and = ^{Ro) is the angular velocity for circular motion at radius 
Rq in the potential Oq. Now consider the situation in the presence 
of the small perturbation d>i. We want orbits that are closed in the 
rotating framej^We expect that far from resonances closed orbits 
will be weak oval deformations of the circular orbits found in Oq. 
We look for closed orbits of the form 

R{t)=Ro^Ri{t), (6) 

0(O = 0oW + 0iW, (7) 

where all quantities with subscript 1 are to be considered small. By 
expressing Eq. ^ in polar components, then substituting Eqs. 


^ Recall that the notion of closureness of an orbit is frame dependent. Or¬ 
bits that are closed in the inertial frame are in general not closed in the 
rotating frame and vice-versa. 


1^ and ([T} into it and approximating to first order in quantities with 
subscript 1, the equations of motion take the following form: 


^1 + 




-q} 


Ri — — 


Ro 


01+2^— 

Rq 


^0 V ^ RoMi) ’ 


V dR 


^0,00 (0 


( 8 ) 

(9) 


where the derivatives of d>i are evaluated along the unperturbed 
trajectory [Rq, %{t)] and are therefore to be considered given func¬ 
tions of time. Exploiting the fact that 0o(O is a linear function of 
time, Eq. ([^ can be integrated immediately to obtain 0i, which can 
then substituted into Eq. i). This gives 


where 




/i(0 = - 


Ri +KqRi — fi (t) , 

(10) 


(11) 

[(''«•■')+ 2a 1 

\dR J Q.fRo ‘ 

(12) 
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Eq.[^is the equation of a forced harmonic oscillator. The natural 
frequency of this oscillator is Kq, the usual epicycle frequency that 
can be calculated from d>o. The driving force is /i(t), and as one 
would expect it is due to the perturbing potential Oi and reduces to 
zero when this is turned off. f\ (t) is a periodic function of time. If 
we expand it in its Eourier components with respect to time, we will 
find the same components contained in the Eourier expansion of Oi 
with respect to 0. This is not surprising, as the perturbed potential 
encountered by the guiding centre of the ballistic particle is given 
by Oi (Rq, 0o(O) where %{t) increases linearly with time. 

The general solution of Eq. {To} could in principle be easily 
written down. Eor each solution, we could eliminate t from Eqs. 

to find the orbit in the form R = R{Q). Not all solutions of 
Eq. |T0] ) correspond in general to closed orbits in the rotating frame. 
Some particular solutions describe closed orbits, while the others in 
general describe non closed loop orbits. The first type of solution 
are the ballistic closed orbits in the epicycle approximation. 

Eor concreteness, let us now consider the particular form of 
the potential used by |Sanders| ( |1977j ) and |Wada| ( |1994j ): 


cFo(R) = 




<&i(i?.e) = 4)bWcos(29), 


(13) 


= -e 


(avoR)'^ 

(if2 + a2)2 ’ 


(14) 


where a, vq, £ are constant parameters. The axisymmetric part d>o 
is a Kuzmin-Toomre potential, and the perturbation Oi was in¬ 
troduced by jSandersj ( |1977j ). Eollowing jWadaj ( |1994j ), we choose 
a = Ikpc and vq = 100km s“^ and a pattern speed of flp = 
10kms“^ kpc“^. Different values of 8 will be considered in the fol¬ 
lowing. Fig. shows the circular speed curve and the behaviour of 
Q. — wKo/ 2, from which we can determine the location of the reso¬ 
nances, for this choice of values of the parameters. For this poten¬ 
tial, Eq. ( p^ becomes 


Ri +KoRi = /wCos(2flft) 


(15) 


where 




V dR 


+ 


2Q. 

Q.fRo 
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(16) 
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Gas flow in barred potentials II 3 



Figure 1. Blue curve: circular speed curve for the potential in Eq. 0 
with values of parameters indicated in the text. Black curves: behaviour 
of — uKq/I for the same potential. The horizontal green line indicates the 
assumed value of the pattern speed. 


The forcing is now a periodic function of time with frequency 2Q.f. 
By solving GD. it is found that ballistic closed orbits have the ma¬ 
jor axis always either perpendicular or parallel to the major axis 
of the bar^The orientation changes abruptly at each ILR, at CR, 
and at OLR. Note however that the above analysis is not valid in 
the vicinity of these points: at all Lindblad resonances the natural 
frequency of the oscillator 2Q.f tends to the forcing frequency Kq, 
while at CR the forcing becomes infinite. The orientation of the 
ballistic closed orbits will be important in the next subsection. 


2.2 Gaseous dosed orbits in the epicycle approximation: 
phenomenological models 


In the previous section we have seen that for the particular potential 
given by Eqs. ( p3] ) and ( p^ the significant ballistic closed orbits are 
always oriented either perpendicularly or parallel to the major axis 
of the bar, with abrupt changes of orientation at each resonance. 
However, the discussion was only valid for ballistic particles. How 
should we modify it to describe the motion of a gaseous parcel 
which is part of a continuous fluid? |Sanders & Huntley|(T976] ) sug¬ 
gested that in the gaseous case the orientation shift between hor¬ 
izontally and vertically aligned orbits happens gradually, and that 
this gives rise to kinematic spiral arms a-la Lindblad. Moreover, 
they noted that the addition of a dissipation term to Eq. HU. 


R\ “h 2^Ri -|- cos 


(17) 


can remove the divergence of the radial oscillation amplitude at the 
Lindblad resonaces, and that the solution of Eq. ( p^ where tran- 


^ These correspond to orbits obtained by taking the limit ^ ^ 0 of the 
solution that neglect transients of Eq. (17} below, which are also the orbits 
with Cl = 0 in the notation of Sect. 3.3.3 in |Binney & TremaineH2QQ8) . 
It should be noted that other closed orbits solutions are possible at special 
radii, i.e. when the ratio Ko/(fl — flp) happens to be a rational number, 
and correspond to orbits with Ci / 0. However these are isolated cases that 
are not part of a family of closed orbits parenting non closed orbits, and 
therefore are not of interest in this paper. 



Rq [kpc] 

Figure 2. Predictions of Wada’s model for the phase (\> as function of 
guiding-centre radius Rq for the potential in Eq. p^ . (]) represents the orien¬ 
tation of the major axis of a closed orbit. Curves corresponding to different 
values of the damping parameter A are shown. The possible values allowed 
by the phenomenological model are enclosed between the blue and yellow 
lines. 


sient^are neglected always describes a closed orbit. They also 
noted that at each ILR, the closed orbit solution so obtained has 
a major axis inclined at 7i/4 independently of the value of X. Since 
7 i/ 4 is halfway between horizontally and vertical elongated orbits, 
this corroborated the conjecture that in the gaseous case we expect 
a gradual shift from vertically elongated to horizontally elongated 
when we cross a Lindblad resonance, which is the key to the gen¬ 
eration of spiral arms in their model. 

|Wada| ( p^^ and [Lindblad & Lindblad] ( p^^ ; |Pinol-Ferrer| 
jet al.j < |20I2j ) independently developed further this idea by taking 
more seriously the dissipation term, considering it as a phenomeno¬ 
logical model that describes the trajectory of gaseous particles over 
extended regions. The approaches of jWadaj ( |I994j ) and [Lindblad] 
[& Lindbl^ ( [1994] ) differ slightly in the way they implement the 
dissipation term. The phenomenological model of [Wada[ ( |l994[ ) is 
based exactly on Eq. GZ}- This amounts to add a dissipation term 
for the motion of the particle in the radial direction only. [Lindblad] 
[& Lindbi^ ( [1994] ) added a dissipation term that is proportional 
to the difference between the velocity of the particle and the local 
circular speed in the underlying axisymmetric potential. Here, we 
consider in more detail the variant of [Wada[ ( |l994[ ), but the results 
of the [Lindblad & Lindblad[ ( p^^ variant are qualitatively similar, 
and a more detailed account of the results using this variant can be 
found in [Pinol-Eerrer et aL] ( [20I2[ ). 

Let us now briefly review Wada’s model. The solution of 
Eq. ITt] ) that excludes transient terms yields closed orbits of the 
following form: 


R(0) =Ro+Acos2(0-4)) , (18) 


^ By transients we mean the part of the solution that vanishes in the limit 
r ^ oo. In this case we have a driven and damped harmonic oscillator. The 
general solution is the sum of a decaying exponential and an oscillatory 
term, and the decaying exponential is the transient. 
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4 Sormani, Binney & Magorrian 



axis. From Eq. we see that the quadrant 2(|) belongs to is de¬ 
termined by the sign of /w, which flips at CR, and by the sign of 
Kq — which flips at each Lindblad resonance. The phase (|) in 
Eq. ( p^ is determined for all orbits once X = X{Rq) is given at each 
radius. Since X has the dimension of a frequency, it is convenient to 
express it as a multiple of a characteristic frequency at that radius. 
Eollowing |Wada| ( |1994| ) , we express ^ as a multiple of the epicycle 
frequency at that radius, namely 

X = Ako . (22) 

Eig. [^reproduces Eig. 2 of |Wada| ( |l994| ) and shows his model pre¬ 
dictions for the phase (|) for different constant values of A. In the 
limit A ^ 0 we recover the ballistic case, while in the limit A ^ oo 
all orbits are circular. The curves A = 0 and A = oo bound the pos¬ 
sible values that the (|) can assume in its model. Note that this plot 
does not depend on the strength of the bar potential 8, because the 
phase (|) does not depend on the magnitude of /w but only on its 
sign. 

A is the amplitude of radial oscillations. Eig. shows the am¬ 
plitude A of the radial oscillations predicted by the model for vari¬ 
ous values of A and two different values of 8. As can be seen from 
this flgure, A always diverges at CR, regardless of the value of A, 
but diverges at the Lindblad resonances only in the absence of the 
damping term. Note also that, for given 8, away from resonances 
the value of A is limited from above: the maximum value that A can 
reach at given Rq is limited by its value for A = 0. Einally, note that 
for a given potential, this theory has only one adjustable parameter 
at each radius, the damping X. Both (|) and A at each radius depend 
on a single phenomenological parameter in this theory. 

Using the theory outlined above, it is possible to construct ex¬ 
plicit models of the spiral arms. FigJS top panels, reproduces the 
model shown in Eig. 4 of |Wada| ( |l994| r It shows a nested sequence 
of closed orbits for values of the parameters 8 = 0.05 and A = 0.05. 
The bottom panels show another model, for a weaker bar 8 = 0.005 
and A = 0.02. 


Figure 3. Amplitude of the radial oscillations predicted by Wada’s model 
for two different values of the bar strength £. The top panel is for £ = 0.05, 
and the bottom panel is for £ = 0.005. In each panel, curves corresponding 
to different values of the damping parameter A are shown. In absence of a 
damping term, A = 0, the amplitude diverges at all Lindblad resonances and 
at corotation. The presence of a damping term makes the amplitude finite at 
the Lindlblad resonances, but not at corotation. 


where the amplitude is 


with 


A = 


|/w| 

F 


(19) 


3 NUMERICAL METHODS 

In our simulations, we assume that the gas is a fluid governed by 
the Euler equations complemented by the equation of state of an 
isothermal ideal gas. Then we run two-dimensional hydrodynam- 
ical simulations in an externally imposed, rigidly rotating barred 
potential. The output of each simulation consists in snapshots of 
the velocity and surface density distributions p(x) and v(x) at cho¬ 
sen times. 

We use a grid-based, Eulerian code based on the second-order 
flux-splitting scheme developed b y |van Albada et ^ 


later used by |Athanassoula| ( p~992| ), [Weiner & Sellwood| ( |19^ I and 


n^i 


and 


others to study gas dynamics in bar potentials. We used the same 


/ r 1 2 


implementation of the code as was used by [Sormani et al. 

(|2015|, 

F=^[K2-(2nf)2] +(4mf)2. 

(20) 

which also implements the recycling law of[Athanassoula[( 

1992). 


We used a grid N xN io simulate a square 8kpc on a side. N 


The phase (|) in (T8) is given by the solution of 
4X\af\ . 


sin 2(|) = 


-sign(/w 


cos2(^ = 


F 


( 21 ) 


that lies in the range — 7U ^ 2(|) ^ 7i. It encodes the information about 
the orientation of the major axis of the closed orbit, and is the incli¬ 
nation of the major axis of the orbit with respect to the horizontal 


depends on the resolution of the simulation. The initial conditions 
of the runs are as follows. We start with gas in equilibrium on cir¬ 
cular orbits in an axisymmetrized potential and, to avoid transients, 
turn on the non-axisymmetric part of the potential gradually. We 
use outflow boundary conditions: gas can freely escape the simu¬ 
lated region, after which it is lost forever. The potential well is suf- 
flciently deep, however, to prevent excessive quantities of material 
to escape. 

The recycling law introduces a term in the continuity equation 
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Figure 4. Top: reproduction of |Wada|jl994) predictions for a bar strength of £ = 0.05. The damping parameter is A = 0.05. Bottom: another model produced 
using |Wad^jl994) theory, for £ = 0.005 and A = 0.02. 


that was originally meant to take into account in a simple way the 
effects of star formation and stellar mass loss. The equation gov¬ 
erning this process is 

^ = a(Po - P^). (23) 

where a = 03 Mq pc“^ Gyr“^ is a constant and po is the initial sur¬ 
face density, which is taken to be po = IMqpc~^. In practice, the 
only effect of the recycling law is to prevent too much gas accumu¬ 
lating in the very centre and to replace gas lost at the boundary due 
to the outflow boundary conditions. It does not affect the morphol¬ 
ogy of the results, so our results do not change if we disable the 
recycling law. 


4 PHENOMENOLOGICAL MODELS VS HYDRO 
SIMULATIONS 

What happens if we run a simulations in the potential of Eqs. HD 
and |T4])? Do the results of a hydro simulation resemble the predic¬ 
tions of the phenomenological models? Fig. [^answers this ques¬ 
tion. It shows the results of a simulation in which E{t) is a slowly 
changing function of time. At ^ = 0 we have 8 = 0, then 8 increases 
linearly with time until it reaches Wada’s value of 8 = 0.05 at 
t = 12.3Gyr, and for later times it is kept constant at this value. The 
potential evolves so slowly that the gas flow conflguration evolves 
almost adiabatically, adjusted to the instantaneous underlying po¬ 
tential. In other words, at each time t the gas conflguration is almost 
the same as the steady state conflguration that would be obtained 
by freezing the potential and then waiting for the gas to settle down 
into a steady state. 

At early times, when 8 is very small, the results of the hy¬ 
dro simulation do qualitatively resemble the predictions of the phe¬ 


nomenological models, as can be seen by comparing the top row 
in Fig. with Fig. At later times, when 8 is greater, this is not 
true: for example panels in the bottom row do not look like |Wad^ 
( |1994| ) predictions. Inspection of streamlines shows that for very 
small values of 8 the gas is flowing on almost circular orbits, as the 
epicycle approximation requires, but when 8 becomes larger the gas 
flows instead on very horizontally elongated orbits, and the epicy¬ 
cle approximation is not valid anymore. How can we explain this 
behaviour? Why does the gas flow on very elongated orbits despite 
the fact that the perturbation is apparently very small? As we shall 
now see, the key to answer these questions lies in the orbital struc¬ 
ture of the underlying potential. 

Fig. shows how the orbital structure of the underlying po¬ 
tential changes as 8 is varied. Each row refers to a particular value 
of 8. Let us first consider the top row, which shows the case 8 = 0, 
when the bar perturbation is turned off and the potential is exactly 
axisymmetric. Consequently, orbits that close in our rotating frame 
are either circular or they are orbits for which the precession fre¬ 
quency happens to coincide with our chosen “pattern speed”, which 
is actually of no physical significance when 8 = 0. The top left panel 
of Fig.j^shows several families of closed orbits: each dot represents 
a closed orbit in terms of its value of the Jacobi constant Ej and the 
coordinate at which the orbit cuts the vertical axis. All the orbits 
shown in this diagram (known as characteristic diagram, see |Con-| 
[topoulos & Grosbol|1989| l have the property of cutting the vertical 
axis at right angles. The line of dark dots shows the circular orbits. 
Since for 8 = 0 the potential is axisymmetric, these are the only 
stable closed orbits. 

The red crosses in the top left panel of Fig. show eccentric 
orbits that close in the rotating frame. Such orbits exist only be¬ 
tween the two ILRs. They spring out from the sequence of nearly 
circular orbits at a radius of the innermost inner Lindblad resonance 
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6 Sormani, Binney & Magorrian 


ILRl, and then merge into it again at a bigger radius corresponding 
to the other inner Lindblad resonance, ILR2. The red crosses that 
are above the dark line are elongated perpendicularly to the bar, 
while those below it are elongated along the bar. In fact, since the 
potential is axisymmetric, the red cross above the dark line and that 
vertically below it represent orbits that differ only in a 90° rotation 
of the major axis. Since the potential is axisymmetric, it is actually 
possible to find equivalent orbits for any orientation, but these are 
not shown in the diagram as they do not cut the vertical axis at right 
angles. The central panel of Fig. shows how the orbits appear in 
the xy plane. These orbits are coloured differently according to their 
value of the Jacobi constant Ej. 

We now consider what happens when a small perturbation is 
turned on, so the potential is no longer exactly axisymmetric. The 
left panel of the second row of Fig. shows that for 8 = 0.005, a 
bifurcation is present approximately at the location of ILRl. The 
line of dark points that marks the circular orbits in the axisymmet¬ 
ric case is now split. The part inside ILRl merged with what used 
to be the orbits elongated along the bar in the axisymmetric case 
(which have become stable) to form a continuous line. In fact, this 
continuous line constitutes the xi family in the notation of |Con-| 
[topoulos & GrosboT] ( |1989| ). The part that used to be circular orbits 
between the two ILRs now forms a closed loop with the orbits elon¬ 
gated perpendicularly to the bar. In fact, the former circular orbits 
between the ILRs have become the stable X 2 family, while the or¬ 
bits elongated perpendicularly to the bar have become the X 3 family 
and have remained unstable, as opposed to the xi orbits that have 
become stable. Note also that the V 3 orbits are not the xi orbits ro¬ 
tated as in the axisymmetric case, because the potential is no longer 
axisymmetric. 

When in the bottom two rows of Fig. we further increase 8 , 
the X 2 -X 2 loop shrinks, and has almost disappeared when 8 = 0.05. 
Note that the shape of the orbits, visible in the central column, does 
change as we increase 8 , but not dramatically. What changes signif¬ 
icantly is the fraction of the volume of phase space that is occupied 
by non-closed orbits that librate around xi orbits rather than around 
X 2 orbits. The right column of Fig. [^illustrates this point by show¬ 
ing surfaces of sectioij^for a value of Jacobi constant Ej = —0.5 
(for a definition of surfaces of section see for example Chapter 3 
in |Binney & Tremain'^|2008| t. When 8 = 0, all non closed orbits 
are parented by the circular orbit, which corresponds to the centre 
of the “eye” in the top-right panel. When 8 = 0.005, two eyes are 
present in the surfaces of section; the centre of the left one corre¬ 
sponds to the xi orbit, the centre of the right one to the X 2 orbit 
(which replaces the circular orbit). Some orbits are now parented 
by the x\ orbit, and some others by the X 2 . As we increase 8 , the x\ 
orbits becomes predominant in the surface of section, and the frac¬ 
tion of orbits parented by it increases until the X 2 orbit disappears 
for this energy. 

It is now easy to go back to Fig. and interpret the results. 
When 8 is small, the gas is circulating on X 2 orbits in the outer 
parts and on weakly elongated x\ orbits in the inner part (inside the 
ILRl)| 3 when 8 is increased, the volume of phase space occupied 

^ However, one must be careful not to confuse the full phase-space volume 
occupied by a group of orbits with the area in a surface of section filled by 
the same orbits, see|Binney et al.|(|1985|. 

^ Note that there must be a transition region in between, and spiral arms are 
generated as the orientation of orbits changes to make this transition. When 
the bar perturbation is stronger and orbits are more elongated, the transition 
can be mediated by shocks (see [Sormani et al.|2015) rather than by a soft 
spiral arm overdensity as in this case. Note also that in the present case the 


by orbits parented by X 2 orbits is reduced. X 2 orbits start disap¬ 
pearing at small radii, and as they disappear the gas has no other 
choice than to settle onto x\ orbits. When 8 = 0.05 the X 2 orbits 
are gone almost everywhere, and the gas is flowing everywhere on 
x\ orbits, including the region where the x\ orbits are highly elon¬ 
gated. In this regime, the epicycle approximation obviously cannot 
work as the x\ orbits are not weakly deformed circular orbits. If we 
call weak a bar that can be well described under the epicycle ap¬ 
proximation, then the 8 = 0.05 case should be classified as a strong 
bar, despite the fact that the non-axisymmetric part of the potential 
might naively appear small compared to the axisymmetric part. 

We have now two questions. 

(i) How do the phenomenological theories reviewed above com¬ 
pare with the hydro simulations in the weak bar case, that is, when 
8 is extremely small and the epicycle approximation is valid? 

(ii) What can we say in the strong bar case, when the gas is 
fiowing on very elongated orbits and the epicycle approximation is 
not valid? 

Sections [^andj^investigate these two questions respectively. 


5 WEAK BAR CASE 

In this section we want to investigate how the phenomenological 
models reviewed above compare in detail with the results of a hydro 
simulation in the weak bar case, when the epicycle approximation 
is valid. To do this, we study the case 8 = 0.005. Note that this is 
a value ten times smaller than the value 8 = 0.05 considered by 
|Wada| < [T9^ and |Pinol-Ferrer et ^ ( |2012| ). As we have seen in the 
last section, that should actually be considered a strong bar case, 
because of its orbital structure. 

When the sound speed is not too high, the morphology of the 
density distribution qualitatively resembles the predictions of |Wada| 
\\99A\ shown in Fig.j^ This can be seen from Fig.j^ which shows 
results of hydro simulations with 8 = 0.005 and different values of 
the sound speed q. In these simulations the bar strength is frozen 
after the value 8 = 0.005 is reached, and the snapshots considered 
are all at t = 2.9Gyr, long after the bar is fully grown to this value 
and the gas has long settled down into a steady-state configuration. 
The left column shows the density distribution. The density distri¬ 
butions for Cs = 0.625 km s“^ and Cs = 1.25kms“^ are very similar: 
it appears as if the gas configuration is tending towards some limit¬ 
ing configuration as q ^ 0 at fixed spatial resolution, and that this 
limiting configuration has much in common with the predictions of 
the phenomenological models. 

A second important prediction of the phenomenological 
model that is confirmed by the hydro simulations is that the spiral 
arms can be understood as kinematic density waves. Fig. [^inves¬ 
tigates the instantaneous streamlines, that is, the streamlines calcu¬ 
lated from a frozen snapshot of the velocity field. Since the gas is 
in a steady state, these are not much different from the real stream¬ 
lines followed by gas parcels during the time evolution. We see that 
in all cases streamlines are almost closed curves, in the sense that 
in general the mismatch after a single loop is much less than the 
extent of a radial oscillation during the loop. 

Another verified prediction is that streamlines are well de¬ 
scribed by closed orbits of the form Eq. ( p~8] ). The middle column in 

transition is X 2 ^ x\ outwards-in, while the transition discussed in |Sormani| 
[et al.||20r5| is x\ -^X 2 . 
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Figure 5. Time sequence for a simulation in the potential defined by Eqs. 03 and 03- Each panel shows the density map at a particular time. The bar strength 
£ increases slowly and linearly with time until it reaches the value of £ = 0.05, and is then frozen. The gas adjusts almost adiabatically to the instantaneous 
underlying potential. The sound speed is q = 2.5kms“^ and the spatial resolution is dv = 20pc. The dotted circles indicate the two ILRs and the CR. In each 
panel, the evolutionary time t and the instantaneous value of £ are shown. The colorbar is in units of Mq pc“^ 
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Figure 6. Orbital structure of the potential given by Eq. (Tg for different values of £. Rows from top to bottom show the cases £ = 0.0,0.005,0.02,0.05, 
respectively. In the left column, the Jacobi constant Ej of a closed orbit is plotted against the value at which the orbit cuts the vertical axis in potentials 
with different values of £. Each dot represents a closed orbit, and all closed orbits cut the vertical axis with purely horizontal velocity. The middle column 
shows the closed orbits in xy space, with orbits of the same colour having the same energy. The orbits shown are those corresponding to values of Ej equal to 
{—0.9, —0.8, —0.7, —0.6, —0.5, —0.4}. The vertical slices of the grid in the left panels are drawn at these energies. The correspondence between energy and 
colours can be identified from the bottom panel in the central column: here greater energy (i.e., less negative) corresponds to bigger orbits, and the same colour 
coding applies to other panels. The right column shows surfaces of section for Ey = — 0.5. 


Fig.j^shows in full lines the same streamlines as the left column in 
the RQ plane, and in dashed lines a simple fit to each streamline us¬ 
ing Eq. ( p^ where A and (|) are considered independent free param¬ 
eters. The streamlines are well fitted by a function of this form, and 
if we were to reproduce a nested sequence of streamlines as in Fig. 
lousing the best fitting values of these two parameters, we would 


obtain kinematic density waves that reproduce those obtained in the 
simulation. 

However, in the phenomenological models, A and (|) both de¬ 
pend on a single parameter, the dissipation X, while in the fitting 
procedure above A and (|) are treated as independent parameters. 
Moreover, in the phenomenological models some values of A and 
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(|) are forbidden: (|) must stay within the blue and yellow line in Fig. 
1^ and the amplitude A away from resonance has an upper bound. 
Therefore it is not obvious that dissipation parameter X required by 
the phenomenological model can be tuned for each streamline as 
to produce the values of A and (|) given by a hydro simulation. In 
other words, the values obtained from a hydro simulation for A and 
(|) need not to come from the same value of X, or could be outside 
the regions obtainable for any value of X. Indeed, from the right 
column of Fig. Hit is clear that in the hydro simulation A and (|) are 
often outside the regions allowed by the phenomenological models 
for the same potential. Therefore not all hydro simulations can be 
reproduced by an appropriate gauging of the dissipation. In fact, we 
chose to fit using A and (|) as two independent parameters instead of 
using directly the dissipation X as the only free parameter because 
A and (|) lie outside the allowed regions too often. However, In the 
limit of very low sound speed, the results of the hydro simulations 
for A and (|) are not far from the predictions of the phenomenolog¬ 
ical model in Fig.|^ The values of (|) in the top two rows of Fig.[^ 
are more similar than those in other rows to the black dots in the 
bottom right panel of Fig. and the values of A have also much 
in common with the blue dots in the same panel, although many 
are still forbidden because they lie just above the curve A = 0 in 
the bottom panel of Fig. Another interesting prediction that is 
verified in this limit is that the orientation of the major axis of the 
closed orbits at the two ILRs is 7i/4. 

Why are the results of the phenomenological models well re¬ 
produced in the limit of vanishing sound speed, and less well when 
the sound speed is higher? To understand this, consider the equa¬ 
tions of motion of a gaseous parcel flowing in the hydro simulation. 
If the fluid were really an ideal isothermal fluid, these would be 


D\ 




(24) 


where D indicates the material derivative, the first term on the right 
hand side is the gravitational force and the second is the pressure 
force. However, in a real hydrodynamical simulation at fixed spa¬ 
tial resolution, some unavoidable amount of numerical viscosity, 
which is not included in Eq. ( [2^ , will always be present. This nu¬ 
merical viscosity decreases as we increase the resolution, and tends 
to zero as the resolution goes to infinity. The equations of motion 
of a parcel of gas are therefore something like 


D\ 


= -V<l.-c2^-Fd^. 


(25) 


where indicates the force arising from the finite resolution, 
which we will loosely refer to as the viscous force. Hence, a parcel 
in a hydro simulation is subject to three different forces: gravita¬ 
tional, pressure and viscosity. 

The phenomenological models completely neglect pressure 
forces. Instead, they only account for the gravitational forces and, 
phenomenologically, for the viscous forces F^. Note that the pres¬ 
sure force can be derived from a “pressure potential” 

<&P = c^log(p/po), (26) 


where po is an arbitrary number. In a steady-state configuration, 
the density does not change with time, and the pressure forces can 
be derived from a static potential. Hence, a gaseous parcel can be 
considered to move in a potential that is the sum of the gravita¬ 
tional plus the pressure potential. This static pressure potential is 
not easily derived a priori, and we can only calculate it once we 
are given the steady-state configuration of the gas after solving the 
hydrodynamical equations. 


The central and right columns in Fig. show the pressure and 
perturbation potential forces along the horizontal axis. The central 
column shows forces in the y direction; the gravitational forces in 
this direction are zero for the chosen form of the perturbation poten¬ 
tial. The right column shows forces in the horizontal direction. We 
see that at high sound speed, when the phenomenological models 
are less accurate, the magnitude of the pressure forces is in general 
greater than the magnitude of the perturbation potential forces. At 
Cs = 1.25 the pressure contribution becomes smaller than the bar 
perturbation contribution, and at q = 0.625 the pressure forces are 
negligible compared to the bar perturbation contribution. This is 
why the phenomenological models above work in this limit. The 
phenomenological models have terms that take into account both 
the effects of viscosity and the gravitational potential, but they do 
not take into account the effects of pressure. In the limit where pres¬ 
sure is negligible, they work. If in the phenomenological models we 
could replace the gravitational potential with the sum of the gravi¬ 
tational plus the pressure one, then they would describe the results 
of the simulations at finite pressure. The problem is that in general 
the pressure potential is not easily obtained a priori. 

Before moving on to discuss the strong bar case, we have a 
couple of remarks left. From a theoretical point of view, when the 
sound speed tends to zero the gas is always supersonic and we ex¬ 
pect all small perturbations to cause shocks. Hence, we expect that 
when Cs tends to zero, we would find shocks anywhere there is a 
small velocity gradient. At a shock, the gradient of the density di¬ 
verges but the forces due to pressure could remain finite as they 
depend on the product of the sound speed (which goes to zero) 
with the gradient of the pressure (which goes to infinity). The rea¬ 
son why shocks are not formed in our simulations when the sound 
speed tends to zero is that at finite resolution the numerical vis¬ 
cosity prevents them from forming. In the limit of vanishing sound 
speed, the numerical viscosity becomes more important than the 
pressure force. Indeed, we have checked that increasing the numer¬ 
ical resolution leads to steady states in which there are sharper den¬ 
sity variations - see Fig.|^ Note that the pressure forces in Fig.j^are 
stronger than the pressure forces for a simulation in the same poten¬ 
tial, with the same value of the sound speed but at lower resolution 
- top row of Fig. |7] Hence, we conclude that the phenomenological 
models work in reproducing the result of a simulation in the limit 
of low sound speed and of finite numerical viscosity, which means 
finite numerical resolution. 


6 STRONG BAR CASE 


In the previous section, we have seen that in the weak bar case the 
spiral arms can be understood as kinematic density waves, gener¬ 
ated by weakly oval closed orbits. What can we say about the strong 
bar case? In order to investigate a realistic example of a strong bar 
we consider a simulation we used in |Sormani et al.| ( [2015| ), per¬ 
formed in the |Binney et al.|(r991| ) barred potential. The simulation 
used here is the one with sound speed q = 10kms“^ and spatial 
resolution dx = 10pc. 

In |Sormani et al.| ( [2015| t, spiral arms were present and noted, 
but not discussed in detail. It was shown that in the region where 
spiral arms are present the velocity field was very well approxi¬ 
mated by x\ orbits. To discuss the spiral arms, we will need to look 
at the tiny deviations of the hydro velocity field from the ballistic 
x\ orbits field. 

The left panel in Fig, shows a density snapshot from the 
simulation in |Sormani et ^ ( |2015| ), taken after the bar is fully 
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Figure 7. Results of a simulation in Wada’s potential in the case of a weak bar (£ = 0.005) for different values of the sound speed c^. Each panel in the left 
column shows the steady-state density long after the bar has reached the final value of £ = 0.005 - compare with Fig.RlThe middle column shows the y 
components of pressure forces (full lines) and bar perturbation gravitational forces (dashed line) from the potential in Eq.]l3|with £ = 0.005 calculated along 
the horizontal axis. The right column shows the x components of the same forces, again calculated along the horizontal axis. Note that pressure forces can be 
derived from a pressure potential <I>p as described in the text. 


grown and the gas has reached an approximate steady state. Spi¬ 
ral arms emerging from the end of the bar are clearly identifiable. 
These spiral arms are stationary in the frame corotating with the 
bar. The central panel of Fig. shows what happens when we 
nest together many instantaneous streamlines. We see that, to a 
very good approximation, streamlines are closed curves that when 
nested together generate the spiral arms. In the right panel ballistic 
closed orbits are shown for comparison. Fig. [TT] analyses in more 
detail a selection of instantaneous streamlines in this snapshot. In 
the top panel, five streamlines are shown in dashed lines. These are 
very nearly closed. In the same panel, the ballistic closed orbits that 
cut the y axis at the same value of the streamlines are shown in full 
lines. It can be seen that streamlines are librations around underly¬ 


ing closed xi orbits. In the middle panel, the same streamlines and 
orbits are shown in the RQ plane. In the bottom panel, the radial 
difference between the two orbits is shown in the RQ plane. The 
librations do not have the simple cos(20) structure that we encoun¬ 
tered in Sect. [5] 

The above considerations strongly suggest that also in the 
strong bar case the spiral arms can be understood as kinematic den¬ 
sity waves, generated by small librations around underlying closed 
orbits. In the epicycle approximation, the librations are around cir¬ 
cular orbits. Here, the gas parcels librate around xi orbits. It is nat¬ 
ural to ask whether the phenomenological models can be extended 
to the strong bar case, where now the perturbations should be made 
around the xi orbits and not around circular orbits. In order to in- 
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Figure 8. Analysis of the instantaneous streamlines for the snapshots shown in Fig.|7] The left column shows a sequence of nested streamlines in the xy plane. 
The middle column shows in full lines the same streamlines in the RQ plane, and in dashed lines the result of a simple fit using Eq. [T^ where A and (j) are 
independent parameters. The right panel shows the best fit values of A and (j) as a function of radius. In each panel the dashed red lines mark the two inner 
Lindblad resonances. 
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Figure 9. Same as the top row of Fig. [7] but for a simulation with twice the spatial resolution. Here, dx = lOpc. 


vestigate this, let us go back to the equations of motion in the ro¬ 
tating frame. For a gaseous parcel in a hydro simulation, Eq. 0 is 
modified to 

X = -V(<& + <&/>) + ilpX - 2Q.p (e^ X x) + Fdj., (27) 

where d>p is the “potential” that gives the pressure force introduced 
above and is the viscous force. We can do a Floquet analysis to 
investigate librations around ballistic closed orbits (see for exam¬ 
ple |Bmney^^^^emame|T987]l. Let us write d> = Op + where 
this time d>o is not necessarily axisymmetric. Let Xc{t) be a ballistic 
closed orbit in the potential Oq. We can write a libration around this 
closed orbit as x{t) = Xc{t) -\-xi{t). By substituting this last equa¬ 
tion into ( [27] ), expanding to first order in quantities with subscript 1 
and cancelling the zeroth order terms, we obtain that the equation 
of motion for the libration is 

Xi =-[(xi-V)V<&o + V$i+V$p]^^(,)+n2xi 

( 2 ») 

-2flp(ej xxO+Fdj • 

In this equation, all the derivatives of the potentials are to be cal¬ 
culated along the unperturbed closed orbit and are therefore known 
functions of time. Equations for the epicycle approximation can 
be considered as a particular case of Eq. ( [28] ), obtained assuming 
that d>o is axisymmetric, Op = 0 and there are no viscosity forces 
(which are then added heuristically in the phenomenological mod¬ 
els). In trying to generalise the analysis of Sect.j^to the case when 
O is a strongly barred potential, we are faced with the difficulty that 
there is in general no natural way of splitting O into a Oq where 
we can easily calculate closed orbits and a perturbation Oi that 
only weakly deforms closed orbits. The only natural choice is to 
set Oi = 0. If this choice is made, the whole perturbing potential is 
given by Op, which is not known a priori. In the phenomenological 
models in the epicycle approximation we were able to neglect Op 
compared to Oi in the low pressure case, but here the same cannot 
be done if Oi =0. 

Suppose that Op is somehow known. If we assume Oi = 0 and 
define a suitable phenomenological viscous term in Eq. ( [28] ), 
would the resulting equation correctly describe librations that give 
rise to spiral arms? To investigate this, we can do as follows. We 
first extract the density from the snapshot shown in the left panel 
of F ig. [To] From this density, we derive the pressure potential using 
Eq. (|26])and Cg = 10 km s“ ^. Then we pretend we do not know how 
Op was obtained and find librations using the following equation: 

Xi = - [(xi • V)V4)o + + n^xi - 2% (Cj X Xi) - 2Xxi 

(29) 

In this last equation, we have introduced a phenomenological dissi¬ 
pation term analogous to the models in the epicycle approximation. 


The dissipation is proportional to the difference between the total 
velocity and the local xi velocity field. d>o, <Fp and are now 
known, and given a value of X we can solve Eq. \29\ to find the 
librations. Among all possible solutions to this equation, we want 
the solutions xi{t) that are periodic with the same period of the un¬ 
derlying closed orbits. This amounts to finding the solution where 
transients are gone, dissipated away. The method that we used to 
solve the differential equation ( |2^ is described in detail in Ap¬ 
pendix 

Fig.|l2|shows a nested sequence of librations around vi orbits 
calculated using Eq. ( |2^ and a value of = 50kms“^ kpc“^ Spiral 
arms identical to those present in the hydro simulation of Fig. {To} 
are recovered. This confirms that the spiral arms can be understood 
as kinematic density waves also in the strong bar case. Eq. ( |2^ well 
describes the librations if the correct Op is provided. 

However, in the above analysis we have cheated in the sense 
that we have obtained Op from the result a full hydrodynamical 
calculation. Is there a simple way of deriving a suitable Op, with¬ 
out solving the full hydrodynaical problem? We haven’t found a 
better way. We considered the possibility of obtaining Op from the 
underlying closed orbits. Since the xi orbits are highly elongated, 
we have tried to infer a density map from orbit crowding, and from 
this the pressure potential. However, the pressure forces obtained 
in this way were much weaker than the one obtained in the hydro 
simulation for a value of the sound speed of Cg = 10kms“^, and 
were not capable of reproducing spiral arms even for a small value 
of the dissipation. For a higher value of the sound speed we ob¬ 
tained segments of spiral arms that were not coherent and did not 
produce any “grand design”. Hence, we conclude that while the 
spiral arms in the strong bar case can still be explained by kine¬ 
matic density waves generated by librations around the appropriate 
closed orbits, calculation of these librations require knowledge of 
the pressure forces and hence solution of the full hydrodynamical 
problem. Indeed, the morphology of the spiral arms found in |Sor-| 
[rnani et ah] ( |20I5| ) depends significantly on the value of pressure, 
and a phenomenological model that wants to reproduce must take 
pressure explicitly into account. 


7 CONCLUSION 

In this paper, we have investigated bar-driven spiral arms in the 
absence of self-gravity. Our focus was on understanding the physi¬ 
cal mechanism involved. We concluded that, both in the weak and 
strong bar cases, the spiral arms can be understood as kinematic 
density waves generated by librations around underlying ballistic 
closed orbits. In the weak bar case, the librations can be considered 
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Figure 10. In the left panel, the density distribution for a simulation from SBM2015 with Cj = lOkrns"' and spatial resolution dx = lOpc. This snapshot is 
taken after the bar is fully grown and the final steady state is reached. The middle panel shows a nested sequence of instantaneous streamlines. The right panel 
shows ballistic closed x\ orbits in the same underlying potential. 


deviations from circular orbits. In the strong bar case, the libra- 
tions are to be considered deviations from the appropriate underly¬ 
ing closed orbits, which can be highly elongated. In the strong bar 
case the epicycle approximation is not valid. In fact, whether a bar 
can be considered weak or strong is determined by the validity of 
the epicycle approximation. In Sect. we argued that a bar is weak 
or strong according as the orbital structure is or is not obtainable 
from the epicycle approximation. Bars that might naively be con¬ 
sidered weak, for example the case 8 = 0.05, should be considered 
strong. 

A parcel in a hydro simulation is subject to three different 
forces: gravitational, pressure and viscosity. We have tested the 
phenomenological models available in the literature aimed at ex¬ 
plaining the spiral arms in the weak bar case, when the epicycle 
approximation is valid. We found that the key ingredient not taken 
into account by these models is pressure. Therefore they work well 
in regimes where pressure can be neglected, such as simulations in 
a weak bar at finite resolution and vanishing sound speed. We have 
also discussed how the phenomenological models should be ex¬ 
tended to the strong bar case. When the pressure forces are known, 
these extensions work very well in explaining the spiral arms in the 
strong bar case. Unfortunately, the pressure forces are in general 
known only after solving the full hydrodynamical problem and are 
not known a priori. Thus, while the phenomenological models pro¬ 
vide insight into the physical mechanism that generates the spiral 
arms, they appear to be of little practical use. 
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panel, streamlines (dashed lines) and ballistic xi orbits that cut the y axis 
at the same value (full lines). In background is visible the density distribu¬ 
tion of Fig.[T^for ease of comparison. The middle panel shows the same 
strea mlin es and orbits in the RQ plane. The bottom panel shows the differ¬ 
ence /?streamiine(0) — ^xi (0)? where the zero is shifted for clarity. 
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APPENDIX A: SOLVING THE FLOQUET EQUATION are periodic with period T\ 


In this Appendix we show how to solve Eq. This equation 
arose during a Floquet analysis to find librations around closed or¬ 
bits in a strongly barred potential. The equation to be solved is: 

Xi = - [(xi • V)V <&0 + + ilpXi - lap (4 xxi)-2Xx. 

(Al) 

Expanding into components, rearranging and renaming we can 
rewrite this equation as the following system: 


Xi = -xiK„{t)-yiKjcy{t)+F^{t)+alxi-2apyi-2kii , 

yi = -xiKxy(t)-yiKyy(t)+Fy{t) + a^pyi+2apXi-2Xyi , 

(A2) 

where 

L 

Fx{t) = — 5 

Fy{‘) = -?y'i>pU)- 

In Eqs. ( |A2| ) it is assumed that > 0 and that all the k’s and F’s are 
given periodic functions of time with period T : 

^xx{t) = Kjcx(^ + T) , 

^xy (0 ~ FT) , 

^yy (0 ~ ^yy FT) , (A4) 

Fflt)=FfltFT), 

Fy(t)=Fy{tFT). 

The goal is to solve Eqs. iA2} to find xi{t) = {xi{t),yi{t)). We 
want to find the solution of this equation that neglects transients, 
i.e., the solution to which all solutions tend in the limit t ^ oo. This 
solution is the one such that xi{t) is periodic with the same peri¬ 
odicity of the F’s and k's. This is also a requirement if we want to 
obtain a closed orbit, and we have to assume such periodicity if our 
phenomenological model is to make sensej^Let us therefore ex¬ 
pand all time dependent quantities as Fourier series, assuming they 


^ Note that Eqs. is conceptually similar to the following simpler equa¬ 
tion 


X + K(r)x -\-'kx = F(t) . (A5) 

If K were constant and not time dependent, Eq. jA5) would be the equation 
of a damped and driven harmonic oscillator. The general solution of this 
equation is a transient that decays exponentially plus a periodic term with 
the same periodicity as F(t). Our problem is more general, and K is a func¬ 
tion of time. We are interested in the solution that neglects the transients, 
which is the generalisation of the solution of the damped and driven har¬ 
monic oscillator that neglects the decaying exponential. This solution is the 
one such that x(t) is periodic with the same periodicity of F(t) and k(t). 


n 

n 

n 

^xy{t) = , (^5^ 

n 

= ^^yy\nd^^^ , 
n 

n 

Fy{t) =ZFy;ne‘"^ ■ 

n 

Fortunately, the product of the two Fourier series gives us another 
Fourier series of the same type. Substituting the Fourier expansions 
into Eqs. ( |A2| ) and then equating coefficients term by term we arrive 
at 

- [(«co)2+n2]x„ = 

Kxx:m^n—m ^.^xv\mXn—m (mCO) |^ 2 AA ;2 + 2 f},pF; 2 j + , 

m m 

-[(«C0)2 + n2]F„ = 

Kyy\rnYn—m y^. ^xv: m^n — m {inCd) |^2X/Lz + Fy-fi . 

m m 

(A7) 

The last two equations are a linear algebraic system of equations in 
the unknowns and F„. All the other Fourier coefficients, F’s and 
F’s, are assumed to be known. To solve this system let us rewrite it 
in the form 


AX = F, 

where A is an infinite matrix and 

/ • \ 


X_1 

Y-i 

Xo 

Yo 

Xi 

Yi 


(A8) 


X = 


: / 


F is 


F = 


(A9) 


(AlO) 


/ : \ 

Fx--i 
Fy--l 
Fx-0 

Fy.O 

FyX 
Fy-l 

\ J 

To write A, let us divide it into a part K containing the Fourier 
coefficients of the k' s and a part P that does not: 


. = P + K. 


(All) 


© 0000 RAS, MNRAS 000, 000-000 






16 Sormani, Binney & Magorrian 


D is a block-diagonal matrix: it has 2x2 blocks along the diagonal. 
Each block is 

^ _ /'-{nay)^-\-2'k{in(o)-\-Q.p {in(o)2Q.p \ 

” —{moy)2Q,p —{n(o)^-\-2X{inQd)-\-Q.pJ 

(A12) 

This is the block referring to the vector {Xn^Yn). The entire matrix 
D is then 




\ 




(A13) 




The matrix K is 




Ko 

Ki 

IK2 

K3 

IK4 





K_i 

Ko 

Ki 

K2 

K3 



IK = 


IK_2 

IK_i 

Ko 

Ki 

IK2 


, (A14) 



K _3 

IK_2 

IK_i 

Ko 

Ki 





K _4 

IK_3 

K_2 

K_i 

Ko 




A- 






'-■) 


where each 

is a 2 

X 2 block given by 







1K„= ( 

K^Kxy,n 

Kxy,n \ 

^yy\n J 



(A15) 


Since the linear algebraic system represented by Eq. ( |A8| l is infi¬ 
nite, we need to truncate it in order to be able to solve it. On the di¬ 
agonal of the matrix K there is Kq. Suppose we truncate the Eourier 
expansion of all k’s at order Kq. Then the matrix A is block diag¬ 
onal and the system is very easy to solve by considering only line 
pairs at a time. Now suppose we truncate the Eourier expansion of 
K.{t) at some higher order. In this case, the matrix A is not block di¬ 
agonal, but is a block band matrix that has elements on the sides of 
the diagonal up to the order of truncation. Eor example, if we trun¬ 
cate at K±i, then we have a (block) band of width 3. To find the 
solution we need to truncate the system at some finite order and in 
general we cannot solve the system exactly as in the case in which 
we truncate at IKq, where equations are decoupled in pairs. Thus 
there are two truncations involved. The first is where to truncate 
the Eourier expansions of K and F. This truncation has to be done 
at some order where K and F are well represented by their Eourier 
expansions. The second is the size of truncation of the system A, 
and must be done after the first truncation has been performed. This 
second truncation must be done for a size sufficiently high that the 
solution is converged and is not affected by a further increase in the 
size of the system. We have found that this convergence takes place 
quite rapidly. 
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